Covid19, Japan
日本の新型コロナウィルスに関する情報が抱える大きな問題は、ソースごとに値が微妙にことなっていることです。原因は明白で、総元締めであるはずの厚生労働省が各自治からの患者データを集計せずに各自治体のプレスリリースやウェブサイトの数値を積上げたものを発表していることにあります。
本来であれば医療機関->市町村->都道府県->厚生労働省の順で患者のデータを全て吸い上げて集計すべきなのです。
Import
厚生労働省
厚生労働省のオープンデータのページから取得した「陽性確定者数」のデータ
df_mhlw <- "https://www.mhlw.go.jp/content/pcr_positive_daily.csv" %>%
readr::read_csv() %>%
dplyr::mutate(date = lubridate::as_date(`日付`),
n = as.integer(`PCR 検査陽性者数(単日)`),
cum = cumsum(n)) %>%
dplyr::select(date, n, cum)
df_mhlw
Covid19Japan
Covid19 Japanが独自に収集している陽性確定者単位のデータ。ソースとデータは全てGitHubにて公開されているが、データはJSON形式である点に注意。発表後に修正されたレコード（インスタンス）は削除されれずにステータスなどが変更されているだけなので、「レコード数 \(\neq\) 累計陽性確定者」である点に注意。
df %>%
dplyr::group_by(patientId) %>%
dplyr::summarise(n = n())df %>%
dplyr::filter(patientId == "-1")df %>%
dplyr::filter(patientStatus == "Deceased")JAG Japan
ジャッグジャパン株式会社が独自に収集している陽性確定者単位のデータ。データはCSV形式にて公開されているが、GIS処理用のデータが含まれている点に注意。「レコード数 \(=\) 累計陽性確定者数」になっている。
df_jag <- "https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv" %>%
readr::read_csv()## Warning: Missing column names filled in: 'X52' [52], 'X53' [53], 'X54' [54]## Warning: 610 parsing failures.
## row col expected actual file
## 4484 市町村内症例番号 a double 神戸59/西宮36 'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## 4485 市町村内症例番号 a double 神戸60/西宮37 'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## 6054 市町村内症例番号 a double 川口40/さいたま50 'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## 12486 市町村内症例番号 a double 神戸237/西宮63 'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## 14454 市町村内症例番号 a double (271)のちに番号なし 'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## ..... ................ ........ ................... ..................................................................
## See problems(...) for more details.df_jag
Area Data
都道府県と地方区分に関するデータをGistで公開していますので、使いたい方はご自由に。
Data Wrangling
ソースによる違い
| 収集組織 | 累計感染者数 | 累計死者数 | 日付 |
|---|---|---|---|
| 厚生労働省 | 98321 | 2020-10-28 | |
| Covid19 Japan | 99674 | 1748 | 2020-10-29 |
| JAG Japan | 93990 | 2020-10-28 | |
| 東洋経済 | |||
| NHK（参考） | 98262 | 1733 | 2020-10-27 |
https://github.com/kaz-ogiwara/covid19/tree/master/data https://www3.nhk.or.jp/news/special/coronavirus/data-all/
Summarize
オリジナルのデータがどのようになっているかskimrパッケージを用いてサマライズしてみます。元がJSON形式なので、読み込んだ直後は殆どの変量（フィーチャー）が文字型になっているのと欠損値が多いことが分かります。
df %>%
skimr::skim()| Name | Piped data |
| Number of rows | 101292 |
| Number of columns | 23 |
| _______________________ | |
| Column type frequency: | |
| character | 19 |
| logical | 3 |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| patientId | 0 | 1.00 | 1 | 8 | 0 | 99675 | 0 |
| dateAnnounced | 0 | 1.00 | 10 | 10 | 0 | 275 | 0 |
| gender | 13323 | 0.87 | 1 | 1 | 0 | 2 | 0 |
| detectedPrefecture | 0 | 1.00 | 3 | 15 | 0 | 49 | 0 |
| patientStatus | 97375 | 0.04 | 8 | 23 | 0 | 8 | 0 |
| notes | 52445 | 0.48 | 1 | 270 | 0 | 46090 | 1 |
| mhlwPatientNumber | 100843 | 0.00 | 1 | 11 | 0 | 434 | 0 |
| prefecturePatientNumber | 11320 | 0.89 | 5 | 20 | 0 | 89963 | 0 |
| prefectureSourceURL | 69996 | 0.31 | 5 | 224 | 0 | 3437 | 0 |
| residence | 21244 | 0.79 | 1 | 38 | 0 | 1421 | 0 |
| sourceURL | 637 | 0.99 | 1 | 239 | 0 | 7804 | 0 |
| relatedPatients | 91195 | 0.10 | 2 | 259 | 0 | 6197 | 0 |
| knownCluster | 98828 | 0.02 | 3 | 88 | 0 | 228 | 0 |
| detectedCityTown | 75911 | 0.25 | 2 | 22 | 0 | 663 | 0 |
| cityPrefectureNumber | 76205 | 0.25 | 1 | 34 | 0 | 25078 | 2 |
| citySourceURL | 89485 | 0.12 | 9 | 317 | 0 | 3627 | 0 |
| deceasedDate | 99534 | 0.02 | 5 | 10 | 0 | 226 | 0 |
| deceasedReportedDate | 100093 | 0.01 | 10 | 62 | 0 | 200 | 0 |
| deathSourceURL | 100232 | 0.01 | 14 | 123 | 0 | 646 | 0 |
Variable type: logical
| skim_variable | n_missing | complete_rate | mean | count |
|---|---|---|---|---|
| confirmedPatient | 0 | 1 | 0.98 | TRU: 99674, FAL: 1618 |
| charterFlightPassenger | 101278 | 0 | 1.00 | TRU: 14 |
| cruisePassengerDisembarked | 101281 | 0 | 1.00 | TRU: 11 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| ageBracket | 0 | 1 | 32.84 | 23.52 | -1 | 20 | 30 | 50 | 100 | ▅▇▅▂▁ |
Transform
各変量（フィーチャー）を適切な形式に変換し、地域区分でも分析できるように都道府県データと結合します。
x <- df %>%
dplyr::mutate(dateAnnounced = lubridate::as_date(dateAnnounced),
ageBracket = forcats::as_factor(ageBracket),
gender = forcats::as_factor(gender),
patientStatus = forcats::as_factor(patientStatus),
residence = forcats::as_factor(residence),
cluster = dplyr::if_else(!is.na(knownCluster), TRUE, FALSE)) %>%
dplyr::select(dateAnnounced, ageBracket, gender, detectedPrefecture,
patientStatus, knownCluster, cluster) %>%
dplyr::left_join(prefs, by = c("detectedPrefecture" = "pref"))
xx %>%
skimr::skim()| Name | Piped data |
| Number of rows | 101292 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| Date | 1 |
| factor | 9 |
| logical | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| detectedPrefecture | 0 | 1.00 | 3 | 15 | 0 | 49 | 0 |
| knownCluster | 98828 | 0.02 | 3 | 88 | 0 | 228 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dateAnnounced | 0 | 1 | 2020-01-15 | 2020-10-29 | 2020-08-11 | 275 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| ageBracket | 0 | 1.00 | FALSE | 13 | 20: 24025, 30: 15110, -1: 13435, 40: 12477 |
| gender | 13323 | 0.87 | FALSE | 2 | M: 49495, F: 38474 |
| patientStatus | 97375 | 0.04 | FALSE | 8 | Dec: 1748, Hos: 1266, Hom: 316, Dis: 283 |
| pcode | 1170 | 0.99 | FALSE | 47 | 13: 31192, 27: 12515, 14: 8698, 23: 6171 |
| 都道府県 | 1170 | 0.99 | FALSE | 47 | 東京都: 31192, 大阪府: 12515, 神奈川: 8698, 愛知県: 6171 |
| 八地方区分 | 1170 | 0.99 | FALSE | 8 | 関東地: 52994, 近畿地: 19858, 九州地: 10997, 中部地: 9831 |
| 広域圏 | 7627 | 0.92 | FALSE | 8 | 首都圏: 53214, 近畿圏: 19297, 中部圏: 8409, 九州圏: 7621 |
| 通俗的区分 | 1170 | 0.99 | FALSE | 11 | 関東: 52994, 関西: 19297, 東海: 8059, 九州: 7621 |
| fct_pref | 1170 | 0.99 | FALSE | 47 | Tok: 31192, Osa: 12515, Kan: 8698, Aic: 6171 |
Variable type: logical
| skim_variable | n_missing | complete_rate | mean | count |
|---|---|---|---|---|
| cluster | 0 | 1 | 0.02 | FAL: 98828, TRU: 2464 |
文字型を因子型に変換するだけでも、変量（フィーチャー）内の水準の比率が大まかにつかめるようになります。
例えば、年代別の陽性判定者数は20代が最も多く、続いて、30代、40代と働き盛りの世代に多いことが分かります。また、都道府県別では、東京、大阪、神奈川、愛知と成っていますが、地方区分別では、以外にも関東、大阪の次に九州がきていることが分かります。
Visualizing
厚生労働省
df_mhlw %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::
ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE)全体傾向
sec_scale <- 50
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced) %>%
dplyr::summarise(n = n()) %>%
tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
to = max(dateAnnounced),
by = "day"),
fill = list(n = 0L)) %>%
dplyr::mutate(cum = cumsum(n),
ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced)) +
ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity",
fill = "dark green", alpha = 0.5, width = 1.0) +
ggplot2::geom_line(ggplot2::aes(y = cum / sec_scale),
colour = "dark green") +
ggplot2::geom_line(ggplot2::aes(y = ma), colour = "dark red") +
ggplot2::scale_y_continuous(
name = "陽性確定数（日別）・7日間移動平均（濃赤折線）",
sec.axis = ggplot2::sec_axis(~ . * sec_scale,
name = "累積陽性確定者数（折線）")
) +
ggplot2::labs(title = paste0("Covid19, Japan（@", Sys.time(), "）"),
subtitle = "", x = "日付")## Warning: Removed 6 row(s) containing missing values (geom_path).
plotly::ggplotly()地方区分別
都道府県別では区分が多すぎるので八地方区分別の積上げ棒グラフを描いてみます。サマリで見たように九州地方が以外にも多いことが読み取れます。
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
to = max(dateAnnounced),
by = "day"),
fill = list(n = 0L)) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) +
ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
ggplot2::scale_fill_brewer(palette = "Set1")
折線グラフで表示してみます。
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
to = max(dateAnnounced),
by = "day"),
fill = list(n = 0L)) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) +
ggplot2::geom_line(ggplot2::aes(y = n)) +
ggplot2::geom_point(ggplot2::aes(y = n), size = 1) +
ggplot2::scale_colour_brewer(palette = "Set1")
さらに分かりやすくするために関東、中部、近畿、九州の四地区を除く他の地区をその他としてまとめてみます。
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
`その他地方` = c("北海道地方",
"東北地方",
"中国地方",
"四国地方"))) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
to = max(dateAnnounced),
by = "day"),
fill = list(n = 0L)) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) +
ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
ggplot2::scale_fill_brewer(palette = "Set1")
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
`その他地方` = c("北海道地方",
"東北地方",
"中国地方",
"四国地方"))) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
to = max(dateAnnounced),
by = "day"),
fill = list(n = 0L)) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) +
ggplot2::geom_line(ggplot2::aes(y = n)) +
ggplot2::scale_colour_brewer(palette = "Set1") +
ggplot2::facet_wrap(~ `八地方区分`, ncol = 2)
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
dplyr::ungroup() %>%
tidyr::pivot_wider(names_from = `八地方区分`, values_from = n) %>%
dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>%
tidyr::pivot_longer(cols = -dateAnnounced,
names_to = "八地方区分", values_to = "n") %>%
dplyr::mutate(`八地方区分` = forcats::fct_relevel(`八地方区分`,
"北海道地方", "東北地方",
"関東地方", "中部地方",
"近畿地方", "中国地方",
"四国地方", "九州地方")) %>%
dplyr::group_by(`八地方区分`) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
fill = `八地方区分`)) +
ggplot2::geom_area(stat = "identity", alpha = 1,
position = ggplot2::position_stack(reverse = FALSE)) +
ggplot2::scale_fill_brewer(palette = "Set3") +
ggplot2::labs(title = paste0("Covid19（@", Sys.time(), "）"),
y = "累計人数")
九州が8月頃から急増しているのは、県別に見ると福岡と沖縄での急増が原因と分かります。
x %>%
dplyr::filter(`八地方区分` == "九州地方") %>%
dplyr::group_by(dateAnnounced, `都道府県`) %>%
dplyr::summarise(n = n()) %>%
dplyr::ungroup() %>%
tidyr::pivot_wider(names_from = `都道府県`, values_from = n) %>%
dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>%
tidyr::pivot_longer(cols = -dateAnnounced,
names_to = "都道府県", values_to = "n") %>%
dplyr::mutate(`都道府県` = forcats::fct_relevel(`都道府県`,
"福岡県", "佐賀県", "長崎県",
"熊本県", "大分県", "宮崎県",
"鹿児島県", "沖縄県")) %>%
dplyr::group_by(`都道府県`) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
fill = `都道府県`)) +
ggplot2::geom_area(stat = "identity", alpha = 1,
position = ggplot2::position_stack(reverse = FALSE)) +
ggplot2::scale_fill_brewer(palette = "Set3") +
ggplot2::labs(title = paste0("Covid19, 九州地方（@", Sys.time(), "）"),
y = "累計人数")
クラスタ比率
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(`八地方区分`, cluster) %>%
dplyr::summarise(n = n()) %>%
tidyr::pivot_wider(names_from = cluster, values_from = n) %>%
dplyr::mutate(ratio = (`TRUE` / (`TRUE` + `FALSE`) * 100) %>% round(1)) %>%
dplyr::rename(`地方` = `八地方区分`,
`非クラスタ感染者` = `FALSE`, `クラスタ感染者` = `TRUE`,
`クラスタ比率[%]` = ratio)Modeling
時系列（TS）分析
日毎の陽性判定者数を時系列データ形式に変換します。判定者数がゼロの日は報告されていないので、最初の陽性判定者が出た日から日日次のカレンダーデータを作成して結合してます。なお、時系列データ形式の周期については日次データなので7日間を周期として設定しておきます（考え方によっては約1ヶ月ごとで4週間=28日を周期にするのもありかと思います）。
tmp <- x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced) %>%
dplyr::summarise(n = n()) %>% print() %>%
dplyr::rename(date = dateAnnounced) %>%
tidyr::complete(date = seq.Date(min(date), max(date), by = "day"),
fill = list(n = 0L)) %>% print() %>%
# tidyr::replace_na(list(n = 0L)) %>% print() %>%
dplyr::select(-date)## # A tibble: 272 x 2
## dateAnnounced n
## <date> <int>
## 1 2020-01-15 1
## 2 2020-01-24 1
## 3 2020-01-25 1
## 4 2020-01-26 1
## 5 2020-01-28 3
## 6 2020-01-29 1
## 7 2020-01-30 4
## 8 2020-01-31 1
## 9 2020-02-01 2
## 10 2020-02-04 1
## # … with 262 more rows
## # A tibble: 289 x 2
## date n
## <date> <int>
## 1 2020-01-15 1
## 2 2020-01-16 0
## 3 2020-01-17 0
## 4 2020-01-18 0
## 5 2020-01-19 0
## 6 2020-01-20 0
## 7 2020-01-21 0
## 8 2020-01-22 0
## 9 2020-01-23 0
## 10 2020-01-24 1
## # … with 279 more rowsts_x <- tmp %>%
ts(frequency = 1)
tsw_x <- tmp %>%
ts(frequency = 7)
tsm_x <- tmp %>%
ts(frequency = 28)# ts_x <- x %>%
# dplyr::filter(!is.na(fct_pref)) %>%
# dplyr::group_by(dateAnnounced) %>%
# dplyr::summarise(n = n()) %>%
# dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>%
# dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>%
# dplyr::select(-date) %>%
# ts(frequency = 1)
#
# tsw_x <- x %>%
# dplyr::filter(!is.na(fct_pref)) %>%
# dplyr::group_by(dateAnnounced) %>%
# dplyr::summarise(n = n()) %>%
# dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>%
# dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>%
# dplyr::select(-date) %>%
# ts(frequency = 7)
#
# tsm_x <- x %>%
# dplyr::filter(!is.na(fct_pref)) %>%
# dplyr::group_by(dateAnnounced) %>%
# dplyr::summarise(n = n()) %>%
# dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>%
# dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>%
# dplyr::select(-date) %>%
# ts(frequency = 28)時系列データに変換したものをプロットすると可視化の項でプロットした棒グラフと同じ形のグラフになることが分かります。
ts_x %>%
plot()
tsw_x %>%
plot()
tsm_x %>%
plot()
上記からトレンド（長期的傾向）を除いたグラフ。デフォルト指定なのでlag = 1
ts_x %>%
base::diff() %>%
plot()
tsw_x %>%
base::diff() %>%
plot()
tsm_x %>%
base::diff() %>%
plot()
トレンド、季節変動（周期変動）、非周期変動に分解した場合。frequency = 1では分解できない点に注意。
tsw_x %>%
stats::decompose() %>%
plot()
tsm_x %>%
stats::decompose() %>%
plot()
トレンドだけを抜き出してみる。
tsw_x %>%
stats::decompose() %>%
.$x %>%
plot(ylim = c(0, 1500))
par(new = TRUE)
tsw_x %>%
stats::decompose() %>%
.$trend %>%
plot(ylim = c(0, 1500), col = "dark green", lwd = 3)
par(new = FALSE)tsm_x %>%
stats::decompose() %>%
.$x %>%
plot(ylim = c(0, 1500))
par(new = TRUE)
tsm_x %>%
stats::decompose() %>%
.$trend %>%
plot(ylim = c(0, 1500), col = "blue", lwd = 3)
par(new = FALSE)
# tsm_x %>%
# stats::decompose() %>%
# .$trend %>%
# plot()スペクトル密度の推定。何を意味するのかよく分からない。
ts_x %>%
spec.pgram()
tsw_x %>%
spec.pgram()
tsm_x %>%
spec.pgram()
自己回帰によるスペクトラム解析
spectrum(ts_x, method="ar")
spectrum(tsw_x, method="ar")
spectrum(tsm_x, method="ar")
Infer
ARIMA Model
ARIMA（Auto Regressive Integrated Moving Average, 自己回帰和分移動平均）モデルによる予測を行ってみます。予測に必要なパラメータはステップワイズにより自動的に最適なものが選択されます。
ts_x %>%
forecast::auto.arima() %>%
forecast::forecast() %>%
plot()
tsw_x %>%
forecast::auto.arima() %>%
forecast::forecast() %>%
plot()
tsm_x %>%
forecast::auto.arima() %>%
forecast::forecast() %>%
plot()